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Deconvolution  in  the  presence  of  additive  noise  is  a  well  known  problem  for  which  there  exists 
a  Wiener  filter  which  simultaneously  spectrally  whitens  while  suppressing  noise.  A  simple  variant  of 
this  standard  Wiener  filter  incorporates  a  parameter,  p  say,  which  is  intended  to  allow  further  weight 
to  be  given  to  noise  suppression.  We  shall  call  such  a  filter  a  modified  Wiener  filter.  To  design  such 
a  filter  it  is  required  to  know  precisely  the  frequency  response  of  the  spread  function  or  wavelet,  plus 
the  spectra  of  the  input  and  additive  noise. 

In  practice  some  response  function  is  taken  to  be  appropriate,  and  the  modified  Wiener  filter 
designed  from  it.  If  the  design  response  function  is  thought  of  as  one  chosen  from  a  set  of  allowable 
response  functions  —  a  realistic  practical  viewpoint  —  then  it  is  shown  how  the  selection  of  the 
design  response,  the  chosen  value  of  the  parameter  p  and  the  noise/input  power  spectral  ratio 
effectively  determine  the  characteristics  of  this  set  of  possible  wavelet  response  functions.  This  is 
demonstrated  for  two  different  error  criteria  —  (i)  the  minimization  of  the  average  mean-squared 
error,  and  (ii)  the  minimization  of  the  maximum  mean-squared  error. 

It  is  shown  how  to  calculate  deconvolution  filters  which  solve  sub-optimal  versions  of  (i)  and 
(ii),  but  which  are  robust  to  uncertainty  in  the  wavelet’s  frequency  response. 
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Deconvolution  by  modified  Wiener  filtering:  interpretation 
for  an  imperfectly  known  wavelet 

A.  T.  Walden 


Introduction 

A  most  important  problem  in  reflection  seismology  is  that  of  deconvolution  plus  filtering. 
Consider  fig.  1 ;  we  can  write 

y(0  =  £  h(t -k)r(k)  +  n(t)  (1) 

k  —  oo 

where  r(t)  is  the  input,  h(t)  is  the  channel  or  spread  function,  n(t)  is  additive  noise 
(uncorrelated  with  the  signal),  and  y(t)  is  the  result.  In  this  paper  we  shall  work  with  the 
discrete- time  representation  (1)  with  unit  sampling  interval.  In  reflection  seismology  we 
would  express  y  in  discrete  time  as  the  familiar  convolutional  model  (1)  written  as 

y(t)  =  h(t)*  r(t)  +  n(t) 

where  we  equate  r(r)  with  the  reflection  sequence,  h(t)  with  the  seismic  wavelet,  and 
n  (/)  with  additive  (colored)  noise. 

We  seek  to  design  a  deconvolution  (or  equalization)  filter  g(t)  which  satisfactorily 
recovers,  or  estimates,  r(t);  see  fig.  1.  There  are  several  statistical  approaches  to  this 
problem,  (for  a  summary  of  other  deconvolution  methods  see  Schultz,  1985).  For  example, 
we  can  carry  out  a  "blind  deconvolution"  where  we  work  only  with  y(t)  and  make  some 
important  assumptions:  (a)  reflectivity  whiteness  and  minimum-phase  wavelet  for  standard 
whitening  or  predictive  deconvolution,  or  (b)  reflectivity  whiteness  and/or  sparse  reflectivity 
for  MED-type  methods.  Both  methods  attempt  to  suppress  the  noise  either  explicitly  or 
implicitly.  For  (a)  a  "prewhitening"  is  applied  to  the  Toeplitz  matrix  before  inversion,  while 
for  (b)  the  harshness  of  the  non-linear  mapping  involved  (see  e.g.,  Walden,  1985)  will  often 
provide  some  protection  against  additive  noise.  A  second  approach  which  is  generally  called 
signature  deconvolution  assumes  that  the  wavelet  is  known  to  a  very  good  approximation. 
This  wavelet  estimate  could  be  model-based,  a  far-field  measured  response,  or  the  result  of 
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other  signal  processing  schemes.  Deregowski  (1971, 1978)  contains  excellent  discussions  of 
signature  deconvolution. 

A  third  approach  involves  Wiener  filtering.  If  the  power  spectrum  of  the  input,  R  (co) , 
and  the  noise,  N  (co) ,  plus  the  frequency  response  of  the  channel  (wavelet)  HT( to)  are  all 
exactly  known,  then  the  mean  squared  error  E  {(r(r)-r(r))2}  is  minimized  when  the 
deconvolution  filter  is  chosen  to  be  the  Wiener  filter 


G(co)  = 


//r(co)/?(co) 

|  Hr(to)  1 2R  (to)  +  N  (co) 


where  HT(t o)  is  the  complex  conjugate  of  HT( co).  Denote  the  class  of  all  possible 
deconvolution-filter  responses  by  & .  Then  G  (co)  is  the  solution  of 


mine  {HT,G) 

& 


where 


e(H,G)  =  E  {(r(t)-f(t))2} 

=  (1/271)/^  {  |  1  -//(co) G (co)  | 2  /?  (co)  +  |G(co)|2N(to)}dco 
Multiplying  top  and  bottom  of  (2)  by  Hj( co) ,  and  then  dividing  through  by  |  HT{ co)  1 2R  (co) 


gives 


G(a»  =  [{ l+p(co)}//r(co)] 


We  have  added  a  superscript  "T"  to  //(co)  to  emphasize  that  //T(co)  is  the  true  response  for 
the  realization  of  the  system  which  is  under  investigation.  Here  p(co)  is  the  noise-to-signal 
ratio  IV  ( co)  /  { |  //r(co)  1 2R  (co)} .  The  amplitude  characteristic  of  the  filter, 

l//r(co)|  R  (co) 

|//r(co)|2/?(co)+N(co) 

provides  noise  rejection  at  the  frequencies  where  R  (co)  is  small  relative  to  N  (co) .  Note  that 
even  a  perfect  all-pass  wavelet  demands  amplitude  as  well  as  phase  compensation. 

Suppose  now  instead  of  merely  minimizing  the  mean  squared  error  E  { (r(r)-r(r  ))2 } 
we  seek  to  minimize 

XjEUrfD-rfr))2}  +  X2E{(n(r)*g(r))2}  (4) 


i.e.,  X]  times  the  mean-squared  estimation  error  +  X2  times  the  mean-squared  filtered 
noise.  The  parameters  X  j  and  X2  are  tradeoff  parameters;  their  relative  sizes  determine 
how  much  effort  is  put  into  minimizing  each  component.  Note  Xi,X2  2  0.  (A  similar 


tradeoff  exists  in  the  Backus-Gilbert  inverse  formalism,  see  e.g.,  Oldenburg,  1984,  p670). 
Such  a  scheme  is  well-known  in  designing  shaping  filters  in  the  presence  of  autocorrelated 
noise,  e.g.,  Robinson,  1980,  p227.  As  shown  in  the  Appendix,  the  modified  Wiener  filter 
resulting  from  the  minization  of  (4)  is 

Gm(g>)  =  [{  1+p  p(o>) } //^(co)]-1  (5) 

where  p  =  (1  +  [X2/Xi])  is  an  adjustable  noise  control  parameter.  Increasing  the  value  of 
p  corresponds  to  increasing  the  desire  for  noise  minimization  at  the  expense  of  spectral 
whitening.  If  X2  =  0,  p  =  1  and  we  recover  (3),  while  if  both  mean-square  components  are 
given  equal  weight,  p  =  2.  For  a  chosen  p  there  are  two  unknowns  in  (5),  viz  p(co) ,  the 
noise-to-signal  ratio,  and  the  frequency  response  of  the  wavelet,  HT{ ©) .  In  this  paper  we 
shall  assume  that  Ht(q&)  is  not  perfectly  known,  but  merely  that  the  estimate  of  the 
wavelet’s  frequency  response,  HD  (o) ,  used  in  the  design  of  the  deconvolution  filter,  has 
been  selected  from  a  set  of  possible  choices.  We  also  assume  that  the  ratio  N  (to)  /  R  (t o)  — 
the  ratio  of  noise  power  spectrum  to  input  power  spectrum  —  is  known.  This  is  not  a 
desirable  assumption,  but  is  necessary  to  confine  all  the  lack  of  information  to  the  wavelet 
In  practice,  using  only  seismic  data,  White  (1984,  pl346)  has  explained,  and  illustrated,  that 
the  noise  and  signal  components  and  hence  the  noise-to-signal  ratio  p(co)  = 
/V (<o) /  {  | HT(a>)  1 (g>)  }  can  be  well  estimated  using  multiple  coherence  analysis. 
(However,  the  estimation  of  HT((o)  from  the  estimate  of  the  signal  |//r(to)| 2 /?(©), 
demands  more  assumptions  of  a  dubious  nature,  e.g.,  white  reflectivity,  and  is  not 
recommended.)  Let  us  denote  our  good  estimate  of  p(to)  by  p(co) .  It  will  be  shown  that  if 
(»)/{  |//D(<o) |2/? (to)}  =  n(co)  say  is  close  to  p(co),  then  results  are  conveniently 
expressible  in  terms  of  p  ,  HD( co)  and  £(co). 

It  will  be  demonstrated  how  the  quantities  p  ,  HD( to)  and  jl(co)  effectively  determine 
the  characteristics  of  the  set  of  possible  wavelet  frequency  response  functions.  This  is 
demonstrated  for  two  different  error  criteria:  (i)  the  minimization  of  the  average  mean- 
squared  error  of  estimation,  and  (ii)  the  minimization  of  the  maximum  mean-squared  error  of 
estimation. 

It  is  also  shown  that  if  HD(to)  is  obtained  by  well-documented  methods,  from  which 
various  error  statistics  can  be  formulated,  then  we  can  design  deconvolution  filters  which 
satisfy  sub-optimal  versions  of  the  two  error  criteria  above,  hence  giving  some  robustness 
against  uncertainty  in  the  wavelet’s  frequency  response. 


Minimizing  average  mean-squared  error  in  theory 

Suppose  the  optimization  problem  is  the  selection  of  the  deconvolution  filter  g(t)  which 
minimizes  the  mean-squared  estimation  error,  averaged  over  all  possible  realizations  of  h(t) 
or  equivalently  the  set  of  all  possible  choices  for  the  wavelet’s  frequency  response,  H( co) . 
Denote  this  class  of  allowable  wavelet  frequency  responses  by  )4r .  Similarly,  denote  the 
class  of  all  possible  deconvolution-filter  frequency  responses  by  &  .  Then  we  seek  to 

min  E^AeiH  ,G)}  (6) 

&  * 

Here  E^.  means  “the  expected  value  over  all  H  in  M'”  and  min  means  “the  minimum 

& 

over  the  class  The  solution  to  this  minimization  problem  is  in  fact  that  given  by  Maurer 
and  Franks  (1970),  namely 


EAH  (co) } 

Gm  =  - *5 - — -  (7) 

£*{  |tf(co)|2}  +  [N (©)/£«»)] 

Comparing  (2)  and  (7)  we  see  that  G  i  in  (7)  is  not  just  the  Wiener  filter  for  the  average 
wavelet,  since  the  denominator  contains  the  term  E^{  \  H  (co)  | 2 }  rather  than 
I £#{//(©)} I2.  Since  Var*{H(c o)}*  £*{ |//(o»|2} -  |£,{tf  (w)}  |2: 

Gi( (n)  =  £^{//*(0))}  /  {  |£*{tf(C0)}|2  +  Var^{H(o»}  +  [ N (©)/£ (co)]}  (8) 

When  the  variation  of  H  is  very  small,  i.e.,  Var^  =  0 ,  the  average  wavelet  is  inverted. 
However,  if  Var^{H{ ©)}  R  (co)  is  comparable  in  size  with  the  noise  term  to) ,  then  the 
characteristics  of  the  optimum  filter  will  be  significantly  influenced  by  the  variation. 

The  modified  Wiener  deconvolution  filter  with  parameter  p  ,  design  response  //D(co) 
and  noise  to  input  power  ratio  N(os)/R  (c o)  has  the  form 

GDm  =  Hoi  co)  /  {  |//D(co)|2  +  p  [N((o)/R{  co)]} 


and  so  has  gain 


=  1  /  {( 1  -i-p  m.(co)  ]  //£,  (co) } 


|Gd(co)|  =  | //£(«) |  /  (|//D(co)|2+p[iV(a»//?(o))]} 


whereas  G  {  has  gain 

I G  i(<0)|=  |£^{//*«o)}|/{  |  £^{//  (£«>)}  1 2  +  Var^{H  (co)}  +  [A/ (CD)//?  (co)] }  (11) 

We  must  now  decide  how  to  relate  HD  (&)) ,  our  designed  response,  to  the  set  of  possible 
outcomes  of  //(co).  It  seems  intuitively  sensible  to  assume  \Hq  (g»I  =  I  £.,{//*  (to)}  I 


because  our  choice  of  Hd(gs)  almost  certainly  involved  "averaging"  using  previous 
experience  of  outcomes,  and  possibly  other  parts  of  the  seismic  section.  So  let  us  take 
I Hq (to) |  —  | E  Then 

Var*  {//(co)}  +  [N(co) //?(©)]  =  p[N(to)/R( 0))] 

i.e., 

Var*{H( to)}  =  (p  - 1)  [N(co)//?(co)] 

=  (p  -l)H(CD)  |//r(C0)|2 

=  (p-l)p(co)  |//D(co)|2  (12) 

=  (p  - 1)  ji(co)  |  Hd  (co)  | 2  provided  p.  =  fi . 

Thus  we  see  that  these  three  quantities  determine  Var  *{H( co)} .  Since  p  =  1  +  [  X2/ X  i  ] 
and  X  j ,  X2  -  0  it  follows  that  p  >  1 ,  and  hence,  as  is  obviously  required,  Var  >  0 .  Note 
that  p  =  1 ,  i.e.,  X 2  =  0 ,  is  equivalent  to  setting  Var*. {H  (co)}  =  0 . 

The  phase  of  the  modified  Wiener  filter  is  -0D(co)  where  Hq(  co) 
=  | Hd (co)  |  e ~ ‘ 0o(co) ,  whereas  the  phase  of  G  i  in  (8)  is  arg  {E *{H* (co)}) .  Thus  we  can 
also  identify 

0d(G»  s  -arg(E*{H\ co)}) 

The  parameter  p  does  not  enter  this  equation;  p  is  only  greater  than  1  when  X2  >  0,  and 
the  term  weighted  by  X2  in  (4)  is  independent  of  phase.  Thus  the  phase  is  unaffected.  Note 
that 

arg  {Hq  (co))  as  arg  {E  *,{H\ CO)}) 

while  we  have  chosen 

|//o(»)l  =  |£*{//*(co)}| 

To  recap,  we  have  shown  that  by  treating  the  wavelet’s  frequency  response  as  a  random 
variable,  and  minimizing  E  *{e(H  ,G)}  over  G  ,  we  get  a  deconvolution  filter  which  is 
identical  to  the  modified  Wiener  filter  designed  from  p,HD( co)  and  N(co)//?(co),  when  we 
set 

|//d(o)I  =  |£*{//*(co)}|;  arg  (H*Dm  =  arg  (E*{H*«x»})  (13a) 

and 

Var  {H (co)}  =  (p  -  l)[N(co)//?(co)] . 

Jr 


(13b) 


Hence  p ,  HD  (co)  and  [N  (to) //?(©)]  can,  in  this  sense,  be  interpreted  as  defining  the 
characteristics  of  the  variation  of  //(co),  while  HD((a )  alone  defines  its  phase.  Remember 
that  the  modified  Wiener  filter  concerns  itself  explicitly  with  the  filtered  noise  (when  p  >  1 ), 
through  the  second  term  in  (4).  The  minimization  of  averaged  mean-squared  error  is  not 
concerned  explicitly  with  filtered  noise,  but  merely  with  mean-squared  estimation  error.  Yet 
the  deconvolution  filters  can  be  made  equivalent  using  (13a)  and  (13b). 


Minimizing  average  mean-squared  error  in  practice 

In  practice  we  won’t  know  E  *,{H*  (ai)}  or  Var^{H( o>)}.  Suppose  instead  we 


obtain  an  estimate  HD( cu)  of  HT( gj)  using  coherence  methods  (White,  1980;  Walden  and 
White,  1984).  Then  we  know 

V  =  Var{HDm  =  £{|//D(0))|2}  -  |£{//D(o))} |2 

=  (l  +  2o2)  \Ht(g»\2  -  |//r(co) |2 

=  2<T2  l//r(CD)  I2 

where  a  is  the  noise  parameter  (see  e.g.,  Walden,  1986  for  more  details).  We  can  estimate 
a2,  and  replace  HT( to)  by  HD( CD)  to  obtain  an  estimate  of  Var  {//^ (co)}  as 

V  =  2a2  |//D(co)|2 

We  now  follow  White  (1984)  to  obtain  an  estimate  jl(to)  of  the  noise-to-signal  ratio,  and 
assuming  |i  =  p  determine  p  from 

1 2 


V  =  (p  -1)  H(CD)  |//D(C0)| 


(i.e.,  equation  (12)  with  Var*,  replaced  by  V ).  Thus  GD(a>)  in  (9)  may  be  computed. 
Instead  of  solving  (6),  GD((0)  will  now  solve 


min  £#.,  {e(H,G)} 
& 


where  M"!  is  the  class  of  wavelet  frequency  responses  with  mean  ///>  (co)  and  variance  V  . 
While  this  is  not  as  ideal  as  solving  (6),  it  is  a  practical  approach,  and  should  prove  to  be 
robust  to  uncertainty  in  the  wavelet’s  frequency  response. 
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Minimizing  maximum  mean-squared  error  in  theory 

In  this  case,  instead  of  the  optimization  (6),  we  consider 


min  max  e  (H,G)  (14) 

6-  * 

i.e.,  first  we  maximize  the  mean-squared  error  over  the  class  of  allowable  wavelet  frequency 
responses,  then  we  minimize  the  result  over  the  class  of  all  possible  deconvolution-filter 
frequency  responses.  The  solution  H' ,  G'  say,  is  a  saddlepoint  to  the  minimax  "game", 
satisfying 

min  e{H',G)  =  e(//',G')  =  ma xe(H,G') 

& 

The  problem  of  determining  the  optimum  G  has  been  studied  by  Moustakides  and  Kassam 
(1985)  who  specified  Pc  to  be  the  class  defined  by 


Al( £0)  £  I // (co)  |  £  ^(co) 


with  0(co)  =  arg  {//(co)}  contained  in  a  known  closed  subset  8  of  (-x,  jc]  (see  fig.  2). 
Moustakides  and  Kassam  determined  that  the  gain,  |  G  2(to)  |  ,  of  the  optimum  filter  has  the 
form 


if  cos{a(co)}>0 
i.e.,  a(o>)£x/2 


-2  cos  {a(to)} 
^z,(w)+i4jy(©) 


I G  2(cd)  l  =  •( 


if  cos  (a(o))}  <  0 
and 

[Al2(o))/?(q»]  2 AL(co) 

Af(co)  “  [Av(o»)-Al(o>)] 


-Al (co)  cos  (a(co)}  if  cos  (a(co)}  <  0 

— 5 -  and 

Al  (g» + [N  (to  )/R  (to))  [4L2(co)  R  (co)]  2  AL  (co) 

N  (co)  <  M  u  (co)  —  (co)] 


(15) 


Suppose  0(co)  is  some  continuous  interval  in  (-x,  x] ,  interpreted  as  a  set  of  points  on  the 
unit  circle  in  the  complex  plane,  then  2a(co)  is  the  angle  subtended  at  the  origin  by  the  arc 
on  the  unit  circle  outside  6(co) ,  see  fig.  2.  If  6(co)  is  a  single  point,  then  a(co)  =  n .  The 
expression  (2x  -  2a(co))/2x  denotes  the  proportion  of  the  circle  corresponding  to  the 
bounds  of  the  uncertain  region  0;  hence  1  -a(co)/x  is  a  measure  of  uncertainty.  Note 
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0(co)  covering  more  that  180°  corresponds  to  cos  (oc(co)}  >  0  and  hence  zero  gain. 

In  (15)  A^(co)  /?  (o>)  is  minimum  signal  power,  so  that  AL2  R/N( co)  is  a  measure  of  the 
minimum  possible  signal-to-noise  ratio  at  co .  Define  a  measure  of  uncertainty  8(oo)  about 
the  channel  gain  characteristic  at  co  by  8(co)  =  [Af/(co)-AL(co)]/[i4y(co)+/l£j(co)] ;  then  the 
degree  of  certainty  =  [l-S(co)]/8(co)  =  2AL (co)/ [A v (co) - AL (co)]  while  the  degree  of 
uncertainty  =  [A£/(co)-AL(co)]/[2AL(co)] .  Thus  we  see  from  (15)  that  if  the  minimum 
possible  S/N  ratio  at  co  is  greater  than  the  degree  of  certainty  of  channel  gain 


I G  2(co)  |  = 


/4L(co)  +  Atf(co) 


i.e.,  the  deconvolution  filter  acts  as  the  inverse  of  the  average  wavelet,  apart  from  the 
attenuation  due  to  phase  uncertainty.  Let  us  compare  (16)  with  the  modified  Wiener  filter 
gain  (10): 

|G0«o)|  =  I  Hq (co)  | /{ | Hd (co)  | 2  +  p  [TV (co)/ R  (co)]} 

In  practice  A^(co)  and  Au{ co)  are  unknown.  However,  if  (say,  90%)  symmetric 
confidence  limits  for  |  HT{ co)  |  have  been  estimated,  they  would  typically  be  of  the  form: 

|//d«b)I(1-G(®))  *  \Ht((0)\  <,  |//D(co)|  (l+(2(co) 

(more  on  this  later).  Equating  AL(co)  with  |  HD  (co)  |  (1  —  Q  (co))  and  /^(co)  with 
|Hd(co)|  (1  +<2  (co))  implies  equating  |//D(co)|  with  -^-(AL(co)  + /^(co) ). 


For  the  gains  to  be  equal,  this  requires: 


-cos  {a(co)> 


\jbm I 

\HDm\2+p  [N (co)//? (co)] 


-cos  (a(co)}  =  1  + 


p  N(  co) 

|  Hd  (co) \zR  (co) 


=  [l+/>n(co)]  .  (17) 

From  (15)  it  is  a  requirement  that  cos{a(co)}  <  0.  Since  [l+pn(co)]-1  >  0  this  will 
always  hold.  Secondly,  |cos{a(co)}|  <,  1  is  a  necessary  condition;  since  p  1 1  and 
(J.(co)  t  0 ,  this  will  be  satisfied.  Suppose  equal  weight  is  given  to  each  component  in  (4), 
i.e.,  Xi=X2,  so  that  p  =  2.  Also  assume  jl(co)  =  jl(co)  =  0.2,  say.  Then 
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cos  (a(co)}  =  -  l/l.4  =>  a(co)  =  71  -  0.775  giving  a  phase  uncertainty  of 
l-2a/2ji  =  1.55/2jc,  which  means  that  @(co)  has  outer  limits  corresponding  to 
approximately  1.55/6.28  of  360°,  i.e.,  about  90°.  The  further  condition  in  (15),  that 
[A*  (os)  R  (co)]/ N  (0))  >  1Al  (as)/  [A  v  (as)  -  AL  (co)]  may  alternatively  be  written 

Av( as)  *  [{2/A?(as)}{N(as)/R(os)}  +  1]  AL(os) 

which  shows  the  role  of  the  noise -to-input  ratio  N(as)/R(as).  Suppose  in  fact  that 
N  (as)/ R  (as)  was  sufficiently  large  that  this  condition  was  not  met,  but  instead  the  minimum 
possible  S/N  ratio  at  co  is  less  than  the  degree  of  certainty  of  channel  gain, 

Av(a»  <  [{2/A?(os)}{N(as)/R(as)}  +  1]  AL(as) 

Then,  providing  cos  (a(co)}  <  0,  the  correct  gain  function  is,  (15), 

-AL(as)  cos  {ot((0)} 


I G  2(co)  |  = 


A?(os)  +  [N(os)/R  (co)] 


(18) 


Comparing  (18)  with  the  modified  Wiener  filter  gain 

|Gd(co)I  =  |//d(co)|/{|//d(co)|2  +  p  [N(as)/R(as)]} 
and  equating  |//O(co)|  with  -^-(A^o^+A^co))  as  before  gives  this  time 


-cos{a(co)}  = 


|ttp(o>i  (A£2(co)  +  [N (as)/R (co)]} 
Al (as)  { I Hp (as) \2  +  p  [N (as)/R (co)]} 


Al(co)  p(co)  I  Hd  (co)  I 


|tfp(co)| 


AL(as) 


[1+P  M-(co)] 


-1 


(19) 


The  second  term  on  the  right  is  identical  to  that  obtained  in  the  case  where  the  minimum 
possible  S/N  ratio  at  co  is  greater  than  the  degree  of  certainty  of  channel  gain,  eqn.(17). 
Clearly  cos  {ot(co)}  is  negative,  since  all  the  terms  on  the  right  in  (19)  are  positive.  What 
about  the  necessary  condition  jcos{a(co)}|  <,  1  ?  By  substituting  |  HD  (co)  | 

=  -^-(A^coJ+Aj^co))  in  (19)  the  condition  is 

|(AL(co)+At/(co))(AL2(co)  +  [/V(co)//?(co)])  S  AL(co){  j(AL(co)+At,(co))2+p  [N (as)/ R (as))} 
which  can  be  written  as 


Al(co)(A^2(co) - [4p  - 2][A/ (as)/R (co)])  £  A{/(co)(AL(co)At/(co)-2[N(co)//?(co)]) 


But  Au(<a)  £  Al  (co)  ,  hence  it  is  sufficient  that 

(Al2(co)  -  [4 p  -  2]  [IV  (to)//?  (to)])  <  (Al  (to)  Av (to)  -2  [N  (to)//?  (to)]) 

Since  Alien)  £  AL{(£>)  A[/((0) ,  it  is  sufficient  that 

[Ap  -2]  [N  (to)//?  (to)]  >  2  [/V  (to)//?  (to)] 
i.e.,  Ap  >  4 
P  >  1 

which  is  always  true. 

The  equation  (18)  shows  that  |  G2(to)  |  acts  as  a  Wiener  filter  for  the  lowest  gain.  An 
additional  attenuation  -cos  (a(to)}  arises  because  of  the  phase  uncertainty. 

The  phase  for  the  modified  Wiener  filter  is  arg  {//^ (to)}  =  -0D(to),  while  the  phase 
for  the  minimax  filter  is,  in  both  cases  discussed  above, 

♦(<*>)  =  n  -  p(to)  (20) 

(Moustakides  and  Kassam,  1985,  eqn.8).  This  phase  is  the  angle  between  0  and  the  line 
OL  which  splits  into  two  halves  the  arc  on  the  circle  outside  ©(to) .  In  order  that  these  two 
expressions  for  phase  be  the  same  we  require 

- 0£> (to)  =  n  -  p(to)  i.e.,  arg  {//o(to)}  =  k  -  p(to)  (21) 

so  that  0o(to)  alone  defines  p(to)  (and  vice  versa)  and  thus  the  position  of  ©(to)  with 
subtended  angle  2a(to) .  This  same  phase  relationship  applies  even  in  the  case  of  a  large 
phase  uncertainty,  when  |  G  2(to)  |  =  0 . 


Minimizing  maximum  mean-squared  error  in  practice 

In  practice  we  won’t  know  the  theoretical  bounds  AL(to),  Ay  (to)  or  the  subset  ©(to)  of 
(-7C,  7t].  However,  we  can  readily  obtain  symmetric  90%  confidence  limits  on  gain  and 
phase  for  the  coherence  methods  mentioned  earlier,  as  AL (to)  =  | HD (to)  |  (l-Q(to)), 
A^(to)  =  |//D(to)|  (l  +  Q(to))  and  ^(to),  ^(to),  (see  e.g.,  Walden,  1986  for  more 
details,  especially  the  form  of  Q  ).  From  fig.  2  the  values  d)L(to)  and  (5*/ (to)  define  an 
estimate  £(to)  of  ©(to) .  If  £(to)  covers  more  than  half  the  circumference  of  the  circle, 
then  the  corresponding  estimate  of  a(to)  satisfies  cos{a(to)}  £  0  and  hence  from  (15)  the 
gain  is  zero.  (The  equivalent  modified  Winer  filter  is  arbitrary,  but  with  p  very  large.) 
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If  cos  (a(co)}  is  found  to  be  less  than  zero,  it  is  necessary  to  test  whether  the  minimum 
possible  S/N  ratio  at  co  is  greater  than  the  degree  of  certainty  of  channel  gain,  i.e.,  we  see  if 

AL2(co)*(co)/A(co)  >  2Al  (to)/ {Au  (co)  -  Al  (co)  } 

But 

AL2(co)/?(co)/iV(co)  =  (l-Q(CO))2  |  ///>  (CO)  I  2  /?  (CD)/  AT  (co) 

=  (1  —  Q  (<o))2/  M-(co) 

As  before,  assuming  |i(co)  has  been  calculated,  and  |i(co)  =  p(co) ,  this  gives 

AL2(C0)/?(C0)/iV(C0)  =  (1  —  Q  (co))2/  tl(co) 

Also, 

2AzXco)/{A„(co-AL(co)}  =  (l-(2(co))/(2(co) 
so  we  must  determine  if 

ji(rn)  £  Q(co)(l-(2(u))  (22) 

Let  us  firstly  take  the  case  where  this  is  true.  From  fig.  2  0L(co)  and  Gy  (to)  define  ($(co) 

A  /A 

and  2a(to) .  Since  0(co)  =  arg  {//0(co)> ,  the  point  estimate  of  the  phase  at  frequency  to, 
is  in  the  middle  of  the  interval  0L(co),  0{/(to) ,  the  condition  (21) 

arg  {Hq  (0))}  =  7t  -  iS(co) 

is  satisfied  automatically. 

By  replacing  p(co)  in  (17)  by  |X(co>  and  a(co)  by  a(co) ,  p  is  defined.  Hence  the 
equivalent  modified  Wiener  filter  is  now  completely  specified.  Instead  of  solving  (14), 
Gd( co)  will  now  solve 

min  max  e(H,  G) 
t- 

where  is  the  class  of  wavelet  frequency  responses  with  bounds  determined  by  the  90% 
confidence  intervals  (A^co),  Au{oS)),  (^(co),  (^(co)).  Again,  while  not  ideal,  this 
approach  should  prove  robust  to  uncertainty  in  the  wavelet’s  frequency  response. 

Of  course  it  is  not  necessary  to  construct  the  modified  Wiener  filter  to  carry  out  this 
minimization,  since  given  A^(co) ,  Ay  (a),  ^  (co)  and  (^(co)  the  minimax  filter  can  be 
constructed  from  (16)  and  (20).  The  modified  Wiener  filter  merely  provides  a  common 
framework.  Note  that  although  the  noise-to-signal  ratio  does  not  enter  (16)  or  (20)  it  is 
needed  to  verify  (22). 


In  the  second  case,  when  }1(03)  >  Q  (to)(l  -  Q  (to)) ,  the  only  change  is  to  find  the  p 
that  solves  (19)  when  a(co),  AL( co)  and  ^I(co)  replace  a(co),  AL( co)  and  |i(co).  As 

before,  |//£>(co)|  =  ^{AL (co)  +  Aa(oi)) . 

Summary 

It  has  been  demonstrated  that  the  modified  Wiener  filter  for  the  deconvolution  problem  can 
be  made  equivalent  to  (a)  the  minimum  average  mean-squared  error  filter,  and  (b)  the 
minimum  maximum  mean-squared  error  filter,  if  certain  quantities  are  equated  in  an  intuitive 
way.  Solving  (a)  and  (b)  requires  perfect  knowledge  of  the  range  of  possible  wavelet 
frequency  responses.  However,  equivalent  sub-optimal  problems  can  be  solved  where 
estimated  uncertainty  characteristics  can  be  calculated  for  the  design  response  used.  The 
deconvolution  filters  for  the  equivalent  sub-optimal  problems  to  (a)  and  (b)  may  all  now  be 
readily  expressed  as  modified  Wiener  filters.  Using  these  methods,  deconvolution  can  be 
made  more  robust  to  wavelet  uncertainties. 
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APPENDIX 

Derivation  of  modification  parameter  p 

Consider  the  minimization  of 

+  X2£{(n(r)*  g(t))2} 

Applying  the  orthogonality  principle: 

*i£{[r(/)-  £  y(t -k)g(k)] -y(j)}  +  X*£{[  £  n{t-k)g(k)]n(j)}  =  0 

k  = i=— • • 

i.e. 

Xl{-Rry(t-j)+  £  Ryy(t-k-j)g(k)}  +  X2{  £  *„„('-* -/)*(*)}  =  0 

k=-~  i  =-■*> 

where  E^fx)  =  £ {r(r)y (r -t)} ;  /^(x)  =  £{r(r)r(r -x)}  and 
£„,(x)  =  £{n(r)n(r-x)}.  It  is  assumed  that  r(r),n(r)  and  hence  y(r)  all  have  mean 
zero. 

Substitute  x  -  t  -j  : 

Xi{-E„(x)  +  £  *„(*"*)*(*)}  +  X2{  £  Rnnfr-VgV)}  =  0 
*=—  *=-- 

Now  Fourier  transform  throughout  to  obtain 

M-S^coHSy/coJCCco)}  +  MS™(g>)G(g>)}  =  0 

where  S^co)  is  the  cross-power  spectrum  between  r(r)  and  y(t),  Syy((o )  is  the  auto¬ 
power  spectrum  of  >(/),  SM((o)  is  the  auto-power  spectrum  of  n(t)  and  G((o)  is  the 
frequency  response  function  of  the  deconvolution  filter  g(t).  Hence 

G (©)  =  X !  S^(co)  /  (X  j  Sw(o»  +  X2  S^f©)} 

Now  y{t)  =  r(/)*  A (r)  +  n(/),  and  the  noise  is  unconelated  with  the  input  Hence 

Sry((0)  =  $,,.(©) //*(©)  =  R((D)H*(G» 

Syy((0)  =  S„(co)  | // (CO)  | 2  +  Sw(©)  =  /?(©)  | //(©)! 2  +  N(o» 

Sn»((0)  =  N(a>) 
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Hence  we  obtain 

G(co)  =  >. ! // ^ (o» /? (to) /  { X j  |//(tO)|2fl(CO)  +  X1N(tO)  +  X2N(<0)} 
=  H  *  (co)  R  (to)  /  { |  H  (to)  | 2  R  (co)  +  [  1  +  (X  2  /X  , )  ]  N  (a» } 

=  H* (a) R (to) /  {\H (fii)\2 R (a)  +  p  N (to)} 
where  p  =  1  +  ( X  2/X, 1 ) 


^*1 


